//
// GPSMapEdit
// (c) Konstantin Galichsky (kg@geopainting.com), 2002-2005
//
// Course over sphere.
//

# include "StdAfx.h"
# include "Course.h"
# include "PtAlgo.h"
# include "Globals.h"

float Course (const point_t & _1, const point_t & _2) {
	if (::fabs (_1.x - _2.x) < 1e-6 && ::fabs (_1.y - _2.y) < 1e-6)
		return 0;

	if (::fabs (_1.y - 90) < 1e-5) // Starting from N pole
		return 180;
	if (::fabs (_1.y + 90) < 1e-5) // Starting from S pole
		return 0;

	const float fX1 = _1.x*c_PI/180.f;
	const float fY1 = _1.y*c_PI/180.f;
	const float fX2 = _2.x*c_PI/180.f;
	const float fY2 = _2.y*c_PI/180.f;

	const float fCosY2 = ::cos (fY2);

	float fCourse = 180.f/c_PI*::atan2 (
		::sin (fX2 - fX1)*fCosY2,
		::cos (fY1)*::sin (fY2) - ::sin (fY1)*fCosY2*::cos (fX2 - fX1)
	);

	if (fCourse < 0)
		fCourse += 360;

	return fCourse;
}

# ifdef _DEBUG
void Test_Course () {
	const point_t pLAX (-(118 + 24/60.f), 33 + 57/60.f);
	const point_t pJFK (-( 73 + 47/60.f), 40 + 38/60.f);

	const float fCourse = Course (pLAX, pJFK);
	assert (::fabs (fCourse - 65.8) < 0.1);

//	const float fCourse1 = Course (0, -90, 0, -1);
}
# endif // _DEBUG